home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / krig2d.pro < prev    next >
Text File  |  1997-07-08  |  10KB  |  285 lines
  1. ; $Id: krig2d.pro,v 1.9 1997/04/28 19:38:43 dave Exp $
  2. ; 
  3. ; Copyright (c) 1993-1997, Research Systems, Inc. All rights reserved. 
  4. ;    Unauthorized reproduction prohibited. 
  5.  
  6. FUNCTION Krig_expon, d, t       ;Return Exponential Covariance Fcn
  7. r = t[2] * exp((-3./t[0]) * d)
  8. z = where(d eq 0.0, count)
  9. if count gt 0 then r[z] = t[1] + t[2]
  10. return, r
  11. end
  12.  
  13. FUNCTION Krig_sphere, d, t      ;Return Spherical Covariance Fcn
  14. r = d/t[0] < 1.0                ;Normalized distance
  15. r = t[2] * (1. - 1.5 * r + 0.5 * r^3)
  16. z = where(d eq 0.0, count)
  17. if count gt 0 then r[z] = t[1] + t[2]
  18. return, r
  19. end
  20.  
  21. FUNCTION krig2d, z, x, y, REGULAR = regular, XGRID=xgrid, $
  22.     XVALUES = xvalues, YGRID = ygrid, YVALUES = yvalues, $
  23.     GS = gs, BOUNDS = bounds, NX = nx0, NY = ny0, EXPONENTIAL = ex, $
  24.     SPHERICAL = sp, NESTED = nest, C0 = C0
  25.  
  26. ;+ 
  27. ; NAME: 
  28. ;    KRIG_2D 
  29. ; 
  30. ; PURPOSE: 
  31. ;    This function interpolates a regularly or irregularly gridded
  32. ;    set of points Z = F(X,Y) using kriging.
  33. ;
  34. ; CATEGORY: 
  35. ;    Interpolation, Surface Fitting 
  36. ;
  37. ; CALLING SEQUENCE: 
  38. ;    Result = KRIG2D(Z [, X, Y]) 
  39. ;
  40. ; INPUTS: 
  41. ;    X, Y, Z:  arrays containing the X, Y, and Z coordinates of the 
  42. ;          data points on the surface. Points need not be 
  43. ;          regularly gridded. For regularly gridded input data, 
  44. ;          X and Y are not used: the grid spacing is specified 
  45. ;          via the XGRID and YGRID (or XVALUES and YVALUES) 
  46. ;          keywords, and Z must be a two dimensional array. 
  47. ;          For irregular grids, all three parameters must be
  48. ;          present and have the same number of elements. 
  49. ;
  50. ; KEYWORD PARAMETERS: 
  51. ;   Model Parameters:
  52. ;    EXPONENTIAL: if set (with parameters [A, C0, C1]), use an exponential
  53. ;             semivariogram model.
  54. ;    SPHERICAL:   if set (with parameters [A, C0, C1]), use a spherical
  55. ;             semivariogram model.
  56. ;
  57. ;   Both models use the following parameters:
  58. ;    A:      the range. At distances beyond A, the semivariogram 
  59. ;          or covariance remains essentialy constant. 
  60. ;          See the definition of the functions below. 
  61. ;    C0:      the "nugget," which provides a discontinuity at the
  62. ;          origin. 
  63. ;    C1:      the covariance value for a zero distance, and the variance
  64. ;          of the random sample Z variable. If only a two element
  65. ;          vector is supplied, C1 is set to the sample variance.
  66. ;          (C0 + C1) = the "sill," which is the variogram value for
  67. ;          very large distances.
  68. ;
  69. ;  Input grid description:
  70. ;    REGULAR:  if set, the Z parameter is a two dimensional array 
  71. ;          of dimensions (N,M), containing measurements over a 
  72. ;          regular grid. If any of XGRID, YGRID, XVALUES, YVALUES 
  73. ;          are specified, REGULAR is implied. REGULAR is also 
  74. ;          implied if there is only one parameter, Z. If REGULAR is 
  75. ;          set, and no grid (_VALUE or _GRID) specifications are 
  76. ;          present, the respective grid is set to (0, 1, 2, ...). 
  77. ;    XGRID:    contains a two element array, [xstart, xspacing], 
  78. ;          defining the input grid in the X direction. Do not
  79. ;          specify both XGRID and XVALUES. 
  80. ;    XVALUES:  if present, XVALUES(i) contains the X location 
  81. ;          of Z(i,j). XVALUES must be dimensioned with N elements. 
  82. ;    YGRID:    contains a two element array, [ystart, yspacing], 
  83. ;          defining the input grid in the Y direction. Do not
  84. ;          specify both YGRID and YVALUES. 
  85. ;    YVALUES:  if present, YVALUES(i) contains the Y location 
  86. ;          of Z(i,j). YVALUES must be dimensioned with N elements. 
  87. ;
  88. ;  Output grid description:
  89. ;    GS:      If present, GS must be a two-element vector [XS, YS],
  90. ;          where XS is the horizontal spacing between grid points
  91. ;          and YS is the vertical spacing. The default is based on
  92. ;          the extents of X and Y. If the grid starts at X value
  93. ;          Xmin and ends at Xmax, then the default horizontal
  94. ;          spacing is (Xmax - Xmin)/(NX-1). YS is computed in the
  95. ;          same way. The default grid size, if neither NX or NY
  96. ;          are specified, is 26 by 26. 
  97. ;    BOUNDS:   If present, BOUNDS must be a four element array containing
  98. ;          the grid limits in X and Y of the output grid:
  99. ;          [Xmin, Ymin, Xmax, Ymax]. If not specified, the grid
  100. ;          limits are set to the extent of X and Y. 
  101. ;    NX:       The output grid size in the X direction. NX need not
  102. ;            be specified if the size can be inferred from GS and
  103. ;          BOUNDS. The default value is 26.
  104. ;    NY:       The output grid size in the Y direction. See NX. 
  105. ; 
  106. ; OUTPUTS: 
  107. ;    This function returns a two dimensional floating point array
  108. ;    containing the interpolated surface, sampled at the grid points.
  109. ;
  110. ; RESTRICTIONS:
  111. ;    The accuracy of this function is limited by the single precision
  112. ;    floating point accuracy of the machine.
  113. ;
  114. ;        SAMPLE EXECUTION TIMES (measured on a Sun IPX)
  115. ;    # of input points    # of output points    Seconds
  116. ;    10            676            1.1
  117. ;    20            676            1.5
  118. ;    40            676            2.6
  119. ;    80            676            7.8
  120. ;    10            1024            1.6
  121. ;    10            4096            5.9
  122. ;    10            16384            23
  123. ;
  124. ; PROCEDURE: 
  125. ;    Ordinary kriging is used to fit the surface described by the
  126. ;    data points X,Y, and Z. See: Isaaks and Srivastava,
  127. ;    "An Introduction to Applied Geostatistics," Oxford University
  128. ;    Press, 1989, Chapter 12.
  129. ;
  130. ;    The parameters of the data model, the range, nugget, and
  131. ;    sill, are highly dependent upon the degree and type of spatial
  132. ;    variation of your data, and should be determined statistically.
  133. ;    Experimentation, or preferrably rigorus analysis, is required.
  134. ;
  135. ;    For N data points, a system of N+1 simultaneous
  136. ;    equations are solved for the coefficients of the 
  137. ;    surface. For any interpolation point, the interpolated value 
  138. ;    is: 
  139. ;          F(x,y) = Sum( w(i) * C(x(i),y(i), x, y)
  140. ;
  141. ;    Formulas used to model the variogram functions:
  142. ;        d(i,j) = distance from point i to point j.
  143. ;        V = variance of samples.
  144. ;        C(i,j) = Covariance of sample i with sample j.
  145. ;               C(x0,y0,x1,y1) = Covariance of point (x0,y0) with (x1,y1).
  146. ;
  147. ;       Exponential covar: C(d) = C1 * EXP(-3*d/A)   if d ne 0.
  148. ;                               = C1 + C0          if d eq 0.
  149. ;
  150. ;       Spherical covar:   C(d) = (1.0 - 1.5 * d/a + 0.5 * (d/a)^3)
  151. ;                               = C1 + C0           if d eq 0.
  152. ;                               = 0                 if d > a.
  153. ;
  154. ; EXAMPLES:
  155. ; Example 1: Irregularly gridded cases 
  156. ;    Make a random set of points that lie on a gaussian: 
  157. ;      n = 15        ;# random points
  158. ;      x = RANDOMU(seed, n) 
  159. ;      y = RANDOMU(seed, n) 
  160. ;      z = exp(-2 * ((x-.5)^2 + (y-.5)^2))    ;The gaussian 
  161. ;
  162. ;     get a 26 by 26 grid over the rectangle bounding x and y: 
  163. ;      e = [ 0.25, 0.0]    ;Range and nugget are 0.25, and 0.
  164. ;                ;(These numbers are dependent upon
  165. ;                ;your data model.)
  166. ;      r = krig2d(z, x, y, EXPON = e)    ;Get the surface. 
  167. ;
  168. ;     Or: get a surface over the unit square, with spacing of 0.05: 
  169. ;      r = krig2d(z, x, y, EXPON=e, GS=[0.05, 0.05], BOUNDS=[0,0,1,1])
  170. ;
  171. ;     Or: get a 10 by 10 surface over the rectangle bounding x and y: 
  172. ;      r = krig2d(z, x, y, EXPON=e, NX=10, NY=10) 
  173. ; 
  174. ; Example 2: Regularly gridded cases 
  175. ;      s = [ 10., 0.2]            ;Range and sill, data dependent.
  176. ;      z = randomu(seed, 5, 6)        ;Make some random data
  177.  
  178. ;    interpolate to a 26 x 26 grid: 
  179. ;      CONTOUR, krig2d(z, /REGULAR, SPHERICAL = s)
  180. ; 
  181. ; MODIFICATION HISTORY: 
  182. ;    DMS, RSI, March, 1993. Written.
  183. ;-
  184.  
  185. on_error, 2
  186.  
  187. s = size(z)        ;Assume 2D
  188. nx = s[1]
  189. ny = s[2]
  190.  
  191. reg = keyword_set(regular) or (n_params() eq 1)
  192.  
  193. if n_elements(xgrid) eq 2 then begin
  194.     x = findgen(nx) * xgrid[1] + xgrid[0]
  195.     reg = 1
  196. endif else if n_elements(xvalues) gt 0 then begin
  197.     if n_elements(xvalues) ne nx then $
  198.         message,'Xvalues must have '+string(nx)+' elements.'
  199.     x = xvalues
  200.     reg = 1
  201. endif
  202.  
  203. if n_elements(ygrid) eq 2 then begin
  204.     y = findgen(ny) * ygrid[1] + ygrid[0]
  205.     reg = 1
  206. endif else if n_elements(yvalues) gt 0 then begin
  207.     if n_elements(yvalues) ne ny then $
  208.         message,'Yvalues must have '+string(ny)+' elements.'
  209.     y = yvalues
  210.     reg = 1
  211. endif
  212.  
  213. if reg then begin
  214.     if s[0] ne 2 then message,'Z array must be 2D for regular grids'
  215.     if n_elements(x) ne nx then x = findgen(nx)
  216.     if n_elements(y) ne ny then y = findgen(ny)
  217.     x = x # replicate(1., ny)    ;Expand to full arrays.
  218.     y = replicate(1.,nx) # y
  219.     endif
  220.  
  221. n = n_elements(x)
  222. if n ne n_elements(y) or n ne n_elements(z) then $
  223.     message,'x, y, and z must have same number of elements.'
  224.  
  225. if keyword_set(ex) then begin       ;Get model params
  226.     t = ex
  227.     fname = 'KRIG_EXPON'
  228. endif else if keyword_set(sp) then begin
  229.     t = sp
  230.     fname = 'KRIG_SPHERE'
  231. endif else MESSAGE,'Either EXPONENTIAL or SPHERICAL model must be selected.'
  232.  
  233. if n_elements(t) eq 2 then begin    ;Default value for variance?
  234.     mz = total(z) / n        ;Mean of z
  235.     var = total((z - mz)^2)/n    ;Variance of Z
  236.     t = [t, var-t[1]]    ;Default value for C1
  237.     endif
  238.  
  239. m = n + 1            ;# of eqns to solve
  240. a = fltarr(m, m)
  241.  
  242. for i=0, n-2 do for j=i,n-1 do begin  ;Only upper diagonal elements
  243.     d = (x[i]-x[j])^2 + (y[i]-y[j])^2  ;Distance squared
  244.     a[i,j] = d & a[j,i] = d             ;symmetric
  245.     endfor
  246.  
  247. a = call_function(fname, sqrt(a), t)        ;Get coefficient matrix
  248. a[n,*] = 1.0            ;Fill edges
  249. a[*,n] = 1.0
  250. a[n,n] = 0.0
  251.  
  252. ; c = invert(a)           ;Solution using inverse
  253. ludcmp, a, indx, even_odd       ;Solution using LU decomposition
  254.  
  255. if n_elements(nx0) le 0 then nx0 = 26    ;Defaults for nx and ny
  256. if n_elements(ny0) le 0 then ny0 = 26
  257.  
  258. xmin = min(x, max = xmax)        ;Make the grid...
  259. ymin = min(y, max = ymax)
  260.  
  261. if n_elements(bounds) lt 4 then bounds = [xmin, ymin, xmax, ymax]
  262.  
  263. if n_elements(gs) lt 2 then $   ;Compute grid spacing from bounds
  264.   gs = [(bounds[2]-bounds[0])/(nx0-1.), (bounds[3]-bounds[1])/(ny0-1.)]
  265.  
  266. nx = ceil((bounds[2]-bounds[0])/gs[0])+1    ;# of elements
  267. ny = ceil((bounds[3]-bounds[1])/gs[1])+1
  268.  
  269. d = fltarr(m)                       ;One extra for lagrange constranint
  270. r = fltarr(nx,ny,/nozero)           ;Result
  271.  
  272. for j=0,ny-1 do begin               ;Each output point
  273.   y0 = bounds[1] + gs[1] * j
  274.   for i=0,nx-1 do begin
  275.     x0 = bounds[0] + gs[0] * i
  276.     d[0] = sqrt((x-x0)^2 + (y-y0)^2) ;distance
  277.     d = call_function(fname, d, t)          ;Get rhs
  278.     d[n] = 1.0                              ;lagrange constr
  279.     lubksb, a, indx, d
  280.     r[i,j] = total(d * z)
  281.   endfor
  282. endfor
  283. return, r
  284. end
  285.